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In this paper we analyze the impact of a circular disc on a free surface using experiments, 
potential flow numerical simulations and theory. We focus our attention both on the study 
of the generation and possible breakup of the splash wave created after the impact and 
on the calculation of the force on the disc. We have experimentally found that drops 
are only ejected from the rim located at the top part of the splash -giving rise to what 
is known as the crown splash- if the impact Weber number exceeds a threshold value 
We crit ~ 140. We explain this threshold by defining a local Bond number Bo t i P based on 
the rim deceleration and its radius of curvature, with which we show using both numerical 
simulations and experiments that a crown splash only occurs when Bo tip > 1, revealing 
that the rim disrupts due to a Raylcigh- Taylor instability. Neglecting the effect of air, 
we show that the flow in the region close to the disc edge possesses a Wcber-numbcr- 
dependent self-similar structure for every Weber number. From this we demonstrate that 
Bo up oc We, explaining both why the transition to crown splash can be characterized in 
terms of the impact Weber number and why this transition occurs for We cr u ~ 140. 

Next, including the effect of air, we have developed a theory which predicts the time- 
varying thickness of the very thin air cushion that is entrapped between the impacting 
solid and the liquid. Our analysis reveals that gas critically affect the velocity of propa- 
gation of the splash wave as well as the time- varying force on the disc, Fry. The existence 
of the air layer also limits the range of times in which the self-similar solution is valid 
and, accordingly, the maximum deceleration experienced by the liquid rim, what sets the 
length scale of the splash drops ejected when We > We cr i t . 



1. Introduction 

The seemingly straightforward experiment of impacting a solid or liquid object on a 
liquid surface exhibits many c hallenges for physical understanding, as was alread y no- 
ticed more than a century ago (| Worthington k Colei[l89l 1 1 90fi IWorthingtoiJl 1 9081) . The 
creation and collapse of an underwater cavity has been studied extensively for a solid im- 
pacting a liquid surface (Richardsonlll948 ^ Gauded[l998l : lLee et aZ.lll997 : Bergmann et al 
l2006UDucleaux et alEoOli iBer gmann et aZ.ll2009i) . as well as the formation of j ets result 



ing from the collapse of this cavity dLonguet-HigginsI 1983 ; Hogrefe et al. Il998t IZeff et al 



2000t iDuchemin et q/j|2002t iGekle et alh00$) . The Trust event that is visible after an ob 



ject hits a liquid is however the splash. It is formed by the liquid that is moving upwards 
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close to the downwards moving object, which initially can be considered as a Wagner 
problem, under the condition th at the object can locally be approximated as flat ( Wagnerl 
19321: IScolan fc Korobkinl 2001). When the object is not locally flat, a splash will only 
be formed provided that certain conditions arc fulfilled. These conditions c an sometime s 
be as subtle as th e wetting properties of a smooth sphere, as was shown by iMav ( 195ll) ; 



Duez et 



( 20071) . 

The impact of a flat plate on a free surface allows for analytical investigations and 
can also be studied using experiments and numerical simulations. S elf-similar solutions 
for th e case without surface tension ( We — > oo) have been found bv llafrati fc Korobkinl 
(l2004h for the i nitial stage after impact, whose scaling exponents were already noticed by 
Yakimovl (|1973l ). Later, this analysis was ex panded in order to calculat e the h ydrodynamic 
load on a flat surface close to impact by llafrati &: Korobkin ( 20081 1 2 1 lh . One of the 
existing differences between our study and the above mentioned ones is that we retain 
surface tension effects in the analysis. We show that self-similar solutions exists for any 
value of the Weber number, a fact corroborated with the aid of boundary integral (BI) 
simulations. We extend our analysis by accounting for air effects, and show that these 
enormously influence both the splash wave velocity during the initial instants after impact 
and the force on the disc. Surprisingly, the gas-to-liquid density ratio also determine the 
sizes of the drops generated when the impact speed is above a certain threshold. Indeed, 
the liquid that is thrown upwards d ue to the im pact of an object, develops in a thin sheet 
with a cylindrical rim on top of it (|Yarinll2006f ) which is susceptible to instabilities that 
may result in the ejection of droplets. Finding the nature of instabilities that result in 
the formation of droplets has been the motivation for many studies, examples of which 
are g iven below. 

In Krechetnikov fc Homsv ( 2009t ). it is argued that the crown formation is the re- 
sult of a Richtmyer- Meshkov instability, and the effect of interfacial curvature on liquid 
rims is discussed in iKrechetnik ov ( 200 91), again i n the light of Richtmyer-Meshkov and 
Ray leigh- Taylor instabilities. Deegan et al. ( 20081) showed that in a narrow range of pa- 
ra meters, this instabil ity can develop in a very regular pattern, which was later used 
by IZhang et~al\ to experimentally show that the wavelength of this regular pat- 

tern agrees with the p redicted value corre sponding to a Raylcigh-Plateau instability. 
In an analytical study, iKrechetnikovl (|2010l ) elaborates more on th e interplay b e tween 
the Rayleigh-Plateau and the Ray leigh- Taylor instability. Recently. iLister et al. ( 2011 ) 
specifically isolated the Rayleigh- Taylor instability on a cylinder, tuning the body force 
normal to its surface by tilting a liquid cylinder inside a liquid of higher density. Finally, 
Lhuissier &: Villermauxl (|2012l ) found that a Rayleigh- Taylor instability is responsible for 
the formation of ligaments from the liquid rim observed in bursting bubbles. 

In this paper we argue that the splash wave deceleration on its top part plays a key 
role in the ejection of droplets. Indeed, we find that drops will be detached when the 
effective weight per unit length of the liquid in the rim, p(A t i P — g) R 2 , overweighs the 
surface tension forces per unit length er, namely, when the local Bond number satisfies 
the condition 



Bo 



tip 



p{A 



tip 



g)R 2 c /o > 1 , 



(1.1) 



where A t i P is the rim deceleration, p and a are, respectively, the density and surface 
tension of water and Rc is the radius of curvature of the rim. We combine experiments, 
numerical simulations and theory to show that the local Bond number at the rim is the 
relevant parameter to predict the transition to crown splash. 

The paper is organized as follows. In Sj2]we experimentally describe the different spatial 
regions appearing after the impact of a solid disc against a free surface and show that the 
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Figure 1. The splash generated at T = 0.01 s after the normal impact against a water interface 
of a circular disc with Rd = 20 mm. (a) A snapshot of the disk just before the impact takes 
place. (6) The splash created by an impact below the critical Weber number (impact velocity 
Vd = 0.6 m/s; We — 99). No drops are ejected from the tip of the rim since We < We cr it- (c) 
The splash created by an impact above the critical Weber number (Vd — 1.0 m/s; We = 274). 
Since We > We cr it, the rim breaks into drops with sizes of the order of the rim width. 



experimental results can be reproduced by using potential flow numerical simulations. 
Making use of potential flow theory, we find in Sj3]that, when air effects are absent, the 
splash region possesses a self similar structure that depends only on We for dimensionless 
times satisfying t <C 1. In section 2] we find the threshold value for the critical Weber 
number We cr a above which the crown breaks up and droplets are ejected. In ij5] we 
reveal the critical role played by air density in setting the maximum deceleration of 
the rim and, using this information, we deduce an expression for the drop size when 
We> We cr i t . In addition, in this section we study the force exerted on the disc during 
the impact process and show how it is regularized by the air layer in between the disk 
and the liquid. Conclusions are drawn in 



2. Experiments and comparison with boundary integral simulations. 

Our experimental setup consists of a disc with radius Rd which we pull down through 
a water surface at a constant speed Vd using a linear motor. The linear motor ensures 
that our disc always moves with a constant prescribed velocity. We record the events 
using a Photron SA1 high speed camera. A m ore detailed description of the experimental 
setup can be found in lBergmann et~al. ( 2009f ). From now on, distances, times, velocities 
and pressures are made dimensionless using Rd, Rd/Vd, Vd and pV£ as characteristic 
length, time, velocity and pressure respectively. Also, variables in capital letters will 
denote dimensional quantities, whereas their corresponding dimensionless counterparts 
will be expressed using lower case letters. 

Figure Q] shows the surface deformation that occurs immediately after a circular disc 
impacts perpendicularly with a constant velocity onto a free surface bounding a deep 
water layer. From these images, it can be appreciated that a circular liquid sheet -the 
splash- is ejected out of the liquid bulk all along the perimeter of the disc. The sheet then 
propagates radially outwards, 'informing' the rest of the fluid of the solid body impact. 
The liquid speeds inside the splash arc much larger than the disc impact velocity Vd- 
indeed, for a given instant in time, the distance traveled out of the liquid bulk by the 
top part of the liquid sheet is much larger than the distance traveled within the liquid by 
the disc. We observe that, depending on the impact velocity, the liquid sheet can either 
breakup into drops or just retract into the liquid bulk without breaking. Indeed, if Vd 
or, equivalently, if the impact Weber number We = pRoV^/a is sufficiently large, the 
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rim at the highest part of the sheet breaks into drops, provoking what we refer to as the 
crown breakup or crown splash [figure QJc)]. If, in contrast, Vd (or We) is sufficiently 
small, surface tension forces and gravity pull back the edge of the rim into the liquid and 
no breakup occurs [figure HJfe)] • 

Note in figures [TH2] that initially the splash wave develops in a confined spatial region, 
which grows in time. More specifically, this region has a characteristic length Xp(T) <C 
Rd located very close to the disc edge, with Xp the distance from the disc edge to the 
point P indicated in figure [5] where the slope in the splash is -1. Figure Eta) shows the 
distance traveled by the point P in the splash wave as a function of time for several values 
of the impact velocity. Clearly, the propagation velocity of the splash wave increases 
as the impact velocity of the disc increases. The different curves in figure 02a), when 
represented using dimensionless units, collapse onto a single one as it is depicted in 
figure meaning that the dimensionless velocity of the splash wave is independent of 
the impact Weber number. The results in figure \%lb) also indicate that the speed of the 
splash wave is initially constant since xp ~ 1.6 — s- Xp ~ 1.6 Vd for t ~ 0. However, after 
a short time interval, the splash velocity is no longer constant but decreases in time. 

Due to the fact that the impact Reynolds number Re = VdRdI v (where v indicates 
the kinematic viscosity of the liquid) is such that Re > O(10 4 ), we expect that viscous 
effects are confined to a thin boundary layer developing at the disc bottom. Thus, the 
velocity field in most of the liquid volume is expected to be described using a velocity 
potential. As a consequence of this, the time evolution of the free surface will be correctly 
predicted using a poten tial flow boundary integral method of the type us ed to describe 
the collapse of cavities ( Longuet- Hiaeins fc Qguall995t iGekle et a/.ll2008r), the ejection 
of Worthington jets (|Gekle et afj l200flt IGekle fc Gordilloll201(A iGordillo fc Gekle l201(]h 
or the formation of bubbles from an underwater nozzle ( Qguz fc Prosperettil 19931 ) . 

Indeed, figure [4] shows a single fluid potential flow simulation of a disc impacting at 
a constant velocity Vd against a free surface which very well reproduces the experiment 
once the origin of times is shifted in time by a quantity to exp . It will be shown below 
that to exp plays the role of a virtual origin in time which accounts for the presence of 
air between the bottom of the plate and the free surface, as will be discussed in 
We conclude from the result depicted in figure 0] that the liquid flow can be accurately 
reproduced for t > to exp using a single fluid potential flow description. 



3. Theoretical description of the splash 

In this section we will provide a theoretical description for the for mation of the splash 
neglec ting the effect of air density, which is in part analogous to llafrati fc Korobkin 
^ct us first notice from figures [T]-[2]that there exist two well-differentiated spa- 
tial regions after a solid impacts a free surface. Indeed, using the (dimensional) polar 
coordinates (R, 8) centered at the disc edge shown in figure [2j we observe that there 
exists an inner region ~ Xp(T) <C -Rd, where the interface deforms appreciably and an 
outer region Xp(T) <C L(T) <C Rd where the interface hardly changes from its initial 
position. We start with deriving the flow field in the outer region, i.e., at distances r from 
the edge of the disc such that xp -C r <C 1, for times close to the moment of impact 
(t <C 1), neglecting air effects and assuming negligible deformations of the free surface 
(Section GO}. After that, we use the analytical expression of this flow field as the far field 
boundary condition for the inner problem: The flow within the splash region, r ~ xp, 
where the deformations of the free interface are no longer small. We will show theoreti- 
cally that this specific type of far field boundary condition leads to the flow field within 
the splash region to be self-similar. We confirm this result by rescaling the splash wave 
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Figure 2. Schematic drawing of the impact of the disc and generation of the splash wave. We 
define the origin of our coordinate system at the edge of the disc for both the cartesian and 
polar coordinate system. The point P, which we use as a reference point to measure velocities 
and lengths in the splash, is defined where the slope of the splash wave is —1. The length of 
the splash region, Xp(T), is the distance between the disc edge and the point P. This variable 
is used to indicate at which distance from the edge of the disc the free surface has deformed 
significantly. We define the intermediate length L(T), such that L(T) ^> Xp(T), but small 
enough when compared to Rd to approximate the flow near the disc edge as 2-dimensional, as 
explained in §[3j 




0.005 0.01 0.015 0.02 0.025 0.03 0.2 0.4 0.6 0.8 1 1.2 1.4 



T(s) t 

Figure 3. (a) Position of the point P defined in figure [2] as a function of time for several values 
of the disc impact velocity, (b) The same curves as in (a) but expressed in dimensionless units. 
Observe that, initially, the speed of the splash wave is constant and proportional to Vd since, 
initially, xp ~ 1.6. After a short time interval, the wave speed is no longer constant but decreases 
in time. 



profiles computed using the numerical boundary integral method for several values of 
the impact Weber number (Section 13.2)) . Finally, we show that the expected self-similar 
scalings related to the vertical positions and vertical velocities are fully recovered after 
we compensate for the downwards motion of the disc, which introduces a non self-similar 
correction (Section [3~3|l . 

3.1. Flow field 

Since xp{t) <C 1, it is possible to define an intermediate length i(t) such that xp(t) <C 
i(t) C 1 at which the height of the interface hardly varies with respect to its initial 
position (see figure [2]). Since l(t) <C 1, the velocity potential 4> at the intermediate region 
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r=2.0 ms 



T=4.0 ms 




7= 6.0 ms 
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Figure 4. Six snapshots of the experiment at different times T after disc impact, compared with 
the corresponding results from the boundary integral simulations (yellow lines). Disc radius is 
20 mm, the impact speed is 1 m/s. The agreement between numerical simulations and exper- 
iments is excellent once two precisions are made: First, in the numerical simulations, the tip 
of the splash is unstable and therefore breaks earlier than in experiments. For this reason, the 
splash in the experiment appears higher than in the simulations. Second, the times correspond- 



ing to the experimental profiles are those of the numerical ones plus a constant Toe 



0.6 ms. 



The existence of this time shift is related to the influence of air at the bottom of the disc and 
can be retained by simply shifting the numerical origin of times by a constant value (virtual 
origin of the self similar solution, see text). 



can be described using a two dimensional approach, which satisfies the following equation: 

2 d 2 4> d 2 cf) 
V ^tfV = °' (3 - 1} 
subject to the following boundary conditions: 

^ = -1 at z = -t~0; x<0, (3.2) 
oz 

which is the kinematic boundary condition imposed by the downward moving disc, 

0-0 at z~0;x>0, (3.3) 
denoting the dynamic boundary condition at the free surfac^J and 

<t> -> for \J X" + 2^ 



oo . 



enforcing the fluid far away from the impact to be at rest. 



(3.4) 



f Equation (|3.3|) has been obtained by approximating the time-integrated unsteady Bernoulli 
equation for very small values of i, which for z = and x > can be written as 
4> ~ 4>(t = 0) - i \V4>(t = 0) | 2 i ~ <t>(t = 0), and taking into account that </>(t = 0) = for 
z = and x > 0. 
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Here, we used the cartesian coordinates (x, z) denned in figure [2] and have taken into 
account the observations in figure H) that for t <§; 1 (T <C 20 ms in figured]), the splash 
develops close to the edge of the disc and that the interface is not appreciably distorted 
in the intermediate region r ~ £(t). The solution to the system (|3.1[) - (|3.4[) can be found 
using standard conformal mapping techniques, yielding the complex flow field 

+ ^.-" 2 «-" 2 (i + (i + ><»r /2 . (3.5) 

where i = y/—l and the constant A is determined by matching the velocity field given in 
equation (|3.5|) with the numerical solution of the tangential velocity field at the disc edge, 
which takes into account the true, i.e., three-dimensional, geometry of the impactor and 
the corresponding flow. Figure [5] shows a comparison for two instants of time, such that 
t <C 1, between the numerical calculated velocity field and the one given by equation (|3.5[) . 
which is valid in the intermediate region r ~ £(t). The agreement is almost perfect when 
the deformation of the free surface is minimal, as in figure [Sf a). Figure shows the 
influence of the splash region r ~ xp, where there is a clear discrepancy between the 
approach of equation (|3.5[) and the numerical solution. However, looking at the region 
r ~ £{t), where deformation can be neglected, the agreement is fully recovered. 

Now, since the solution (|3.5[) constitutes the outer boundary condition for the velocity 
field in the inner region r ~ xp (t) <C £(t) < 1, we are interested in the approximation of 
this equation in the limit r«l, which yields, 

l^-Ar-^sin^), ^--l-Ar-^cosieM (3.6) 
and the potential 

</>~-z-2Ar 1/2 sin(6»/2), (3.7) 

which will serve as a boundary condition to the problem of determining the free surface 
shape and potential in the splash region r ~ xp(t). Notice from the first equation in 
(|3.6[) that the modulus of the velocity component tangent to the disc bottom surface 
(8 = 7r) diverge as r -1 / 2 . We will show that the huge tangential velocities near the disc 
edge, which have to accommodate to the mass of fluid at rest, are responsible for the 
generation of the splash wave. 

3.2. Self similarity 

In the frame of reference moving at the disc velocity, the splash region r ~ Xp(t) <C 
£(t) can be described by solving the Laplace equation fl3.1[) subjected to the following 
boundary conditions at a given instant in time: 

d<t> , s 

--=0 at z = x < , 3.8) 
oz 

which is the kinematic boundary condition at the bottom of the disc, and 

which is the dynamic boundary condition at the free surface, with We = pV^Rd/o the 
Weber number and Fr = V^/lgRo) the Froude number. These equations need to be 
complemented with the far field velocity potential in the region where the interface is 
virtually undisturbed, i.e., 

4> -> -2Ar 1/2 sin(0/2) for r -> oo . (3.10) 
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Figure 5. The analytical flow field given in equation (|3.5[) compares very favorably with the 
numerical one once the constant A is set to 0.44. (a) At an extremely short time after impact, 
when the free surface has not yet deformed, the analytical solution agrees with the boundary 
integral result in the full inner domain where r < 1. (6) At a later point in time we observe that 
close to the splash region, where the deformation is appreciable, the analytical flow field deviates 
from the numerical solution. Away from the splash region, the agreement improves again. 

The function f(x,t) defines the position of the free interface, which satisfies the kinematic 
boundary condition 

dj_ = di _ d^di 

dt dz dx dx 

hlwhere 

f{x, t = 0) = and f(x -> oo, t) -> t. 
Note that in the Bernoulli equation (|3.9[) we have used the interfacial curvature 



at z = f(x, t) x > , 



(3.11) 
(3.12) 
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to express the pressure jump across the surface, i.e., p = nWe~ . Since there is no 
characteristic length scale in the system of equations (|3.1|) and (|3.8|) - f|3.11|) . we expect 
the existence of self-similar solutions of the type 



<f> = (t - t ft 



which we now write as 



(t-to) a, (t-t ) a 

M x,v), 



(3.14) 



t"<P{X,V), (3-15) 
t — tg. The dimcnsionlcss quantity to can be an 



with x — X / Ta > V — Z /T a y an d T 
arbitrary constant, and the shape of the free surface can be expressed as 

f(x,t)=T a F( X ). 



(3.16) 

From the Bernoulli equation (|3.9|) we find, by comparing the first two terms, that self- 
similar solutions can only exist if j3 = 2a— 1. Moreover, the matching with the boundary 
condition at infinity (|3.10p imposes the additional restriction 2/3 = a. Combining these 
two conditions then results in j3 = 1/3 and a = 2/3. Thus, lengths are expected to scale 
with t 2 / 3 and velocities with r -1 / 3 . Indeed, the system of equations that solve for both 
4> and F reads, with relative errors ~ (^(t 1 / 3 ) <C 1, 



d 2 4> 
dx 2 

on 



d 2 4> 
drj 2 



0, 



at n = , x < 



-2Af 1 / 2 sin(6»/2) for 



with 



at n 



d<j) , dej) 
dr) 



ox 
F(X), 



do 



do 

dr] 



K 

~We 



vV + X 2 

r 4/3 



Fr 



= 



(3.17) 

(3.18) 
(3.19) 

(3.20) 



with 




-3/2 



and 
2 



dF 



-F - x- 

3 dx 



dr] 



dFdcp 
d X dx 



at n = F(x) and F(x) 



r2/3 



for 



oo . 



(3.21) 

Note that there are two terms breaking the self-similarity of equations (|3.17[) - (|3.21[) 
namely, the last term of equation (|3.20[) and the far field boundary condition in equation 
(|3.21[) due to the presence of t in them. There will therefore exist a self-similar solution 
whenever r > t and both r 1 / 3 <C 1 and r 1 / 3 <C Fr 1 ^ 4 (~ 0(1)) are satisfied. In other 
words, assuming that to is sufficiently small, the self similar solution will be valid during 
a short time interval after impact. Under the conditions in which the flow field is self 
similar, the disc has moved down only slightly, the normal velocity boundary condition 
is then approximately applied at the undisturbed free surface and gravitational effects 
can be neglected. 

It is however crucial to notice that the Laplace pressure term nWe^ 1 in equation 
(|3.20p is self-similar. This is because the result of balancing the Laplace pressure term 
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Figure 6. (a) Time evolution of the horizontal and vertical coordinates of point P. Numerical re- 
sults reveal that the scaling for both x v and z p follows the prediction of equation equation (|3.14|l 
for a = 2/3. As explained in the main text, the results for z p slightly deviate from a pure power 
law due to the real boundary condition at infinity in (|3,10[) . (&) Comparison between the exper- 
iment and the simulation for the radial position of point P. The green diamonds indicate the 
unmodified experimental data. Once a time shift i.e, a virtual origin for times to = toexp — 0.03 
is introduced, numerics (blue circles) and experiments (red dots) are in excellent agreement. 
Deviations from the scaling are expected when t ~ 1, because of geometrical and gravitational 
effects that break the self-similar scaling. 



in the Bernoulli equation (|3.20j) with the inertial terms on the left hand side using the 
self-similar ansatz (|3.15p - (|3.16p will give us 2(3 = a and (3 = 1 — a, which is solved by 
exactly the same exponents a = 2/3, /3 = 1/3 that we have just found by matching 
to the far field boundary condition (|3.10[) . It is this remarkable feature that warrants 
the existence of self-similar solutions to the system (|3.17[) - (|3.2ip for every value of the 
Weber number We. To check our theoretical predictions, figure |6] represents the function 
xp(t) calculated using the single fluid potential flow numerical code finding that, as 
expected, xp(t) oc f 2 / 3 for t<l. For times t ~ 0(1) the t 2 ^ 3 self-similar behavior breaks 
because the non self-similar terms in equations (|3.17p - (|3.21[) are no longer negligible. The 
experimental data included in figure |6j&) shows that the radial position of point P can be 
fitted with a function proportional to (t — toexp) 2 ^ 3 , with to exp ~ 0.03 the virtual origin 
of times, which accounts for physical effects not taken into account in the theoretical 
derivation and the numerical simulations. 

Figurc[BJa) shows that, contrary to xp(t), the time evolution of the vertical coordinate 
of point P calculated numerically, zp(t), is not purely proportional to t 2 / 3 . As will be 
shown in § 13. 3[ this is due to the non self-similar effect of the downward motion of the 
disc or, equivalcntly, the upward motion of the undisturbed free surface in the frame of 
reference moving with the disc. 

An alternative way of understanding the xp oc r 2 / 3 scaling is obtained by making use 
of the fact that the boundary condition at infinity (|3.19p implies that the disc acts as a 
constant 'source' of momentum in time. Indeed, making use of the velocity field calcu- 
lated from the potential (|3.19p . the variation per unit time of the horizontal momentum 
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Figure 7. The surfaces E_r (semicircle of radius L <C Xp), Ed (line of length L located below 
the disc) and the free interface are the boundaries limiting the volume fl(f) to which the integral 
balance of horizontal momentum is applied. 

enclosed in the volume fi denned in figure ([7]) satisfies the following integral balance: 

e x ■ — [ pvduj = e x ■ [ [ pvv • nR d9 + p Vd [ v dx — / pnR d6 | = 
dtJn \J n Jo h J 

f 2rr f Ro 86 

-pA 2 VlR D \ sin(0/2) (sin(0/2)cos0 - cos(0/2)sin0) d9 + pV D / -jr- dx+ 

Jtt JO U x 

1 C 27r 7T f R ° 86 TT 

+ - P A 2 V%R D cos6d6 = p-A 2 Vt > R D +pV D / dx = p-A 2 Vg R D + 

+ O(pV£r d (R /Rd) 1/2 ), 



o 



(3.22) 

where n is the outward normal vector to the surface defined in figure El (e x ,e z ) are 
used to indicate the unit vectors along the coordinates (x, z) and use of the steady Euler- 
Bernoulli equation namely, p+l/^pv-v = C, has been made to express the pressure 
field p over the surface of figure [7] as a function of the velocity field. Notice that the 
contribution of the dimensionless flux of momentum through the surface defined in 
figure [7] is of the order of (Rq/Rd) 1 ^ 2 *C 1 and, therefore, can be neglected during the 
initial instants of the splash. 

The result in equation (|3.22p expresses that the time variation of the horizontal mo- 
mentum contained in the volume Cl(t) is constant. We could have used this argument to 
deduce that there exists a self-similar solution for short times after impact. Indeed, the 
characteristic length of the spatial region where the flow is unsteady is ~ Xp(T) and, 
consequently, velocities within the splash region should scale as V oc Vd xp/r. Therefore, 

d f ^ 
i • — / pvduj oc pVlRn^ ^pVl Rd^A 2 -> x P oc r 2 / 3 , (3.23) 
at J n t 2 2 

which concludes our alternative derivation of the scaling law for xp. 

Clearly, the solution of the system (|3.17j) - (|3.21j) needs to be found numerically. It 
is, however, easier to solve the Laplace equation subjected to the unsteady boundary 
conditions given by (|3.9|l - (|3.11|> and then express the solution in terms of the variables 
X, rj and F defined in equations (|3.15[) - (|3.16[) . i.e., by just using the boundary integral 
method. 

Figure [S] shows the solution from the boundary integral method for one specific Weber 
number at three instances in time, which all are in the regime r <C 1 where we expect 
to find self similar solutions. In figure |9jc) we have rescaled the solution of figure [U 
according to (|3.15j) . Inspecting all shapes in figure [9^ a- d), we see that indeed we find self 
similar solutions for a large range of Weber numbers, that only depend on We. Figure ITTJl 



12 



Ivo R. Peters, Devaraj van der Meer and Jose Manuel Gordillo 



0.5 



We =21 A 




x 



Figure 8. Time evolution of the splash region for a given value of the Weber number. The 
origin of the coordinate z is located at the disc position. Notice that the height of the splash 
wave is much larger than the depth at which the disc is located. Considering a disc with a radius 
Rd = 20 mm and impacting at a speed Vd = 1 m /s, the dimensionless times shown in this 
figure, t = 5 x 10~ 3 , r = 4 x 10~ 2 and r = 10 _1 , correspond to T = 0.1 ms, T = 0.8 ms and 
T — 2 ms respectively. 



sh ows that for We — > oo the solution becomes independent of We, confirming the results 
of llafrati k Korobkinl <|2004h . 



3.3. Correction for non self-similar terms 

From figure [6] it is clear that the vertical position z p of point P does not follow the scaling 
r 2/3 gtrictly^ except in the limit r — > 0. This difference is caused by the downward motion 
of the disc (or, equivalently, the upward motion of the undisturbed free surface, since 
the problem (j3.17B3.2TI) was formulated in the frame of reference comoving with the 
disc), which introduces a non self-similar term in the far field boundary condition of 
equation (|3.21[) . Nevertheless, we can compensate for the effect of the downward motion 
of the disc by introducing a constant B z such that 

(z p + B z t)(xt 2 ^. (3.24) 

A similar correction is needed for the vertical velocity u p at point P: 

(u p + B u ) oc r" 1 / 3 . (3.25) 

Figure QTJ shows the position, velocity, and width of the splash at point P. All values 
become independent of time after the proper rescaling and figure [T^] shows the scaling 
of lengths (a) and velocities (b) for a wide range of Weber numbers. As expected from 
the previous analysis, the rescaled values become independent of We for large Weber 
numbers. The same holds for the correction terms: B z = 0.56±0.03 and B u = 0.59±0.01 
defined in equations (f3~24l) - ([3~25l) for We > 50. 



4. Crown breakup transition 

We have provided a scaling for the shape of the splash, and have shown that the splash 
possesses a self-similar shape for every value of the Weber number. We will now have a 
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Figure 9. Shape of the splash with all distances rescaled by r 2 ^ 3 for four different values of the 
Weber number We. Each plot contains three instances in time (r = 0.005, r = 0.04, r = 0.1). 
The shapes have first been translated such that point P is in the origin of the plots, in order to 
better show the collapse of the different shapes. 
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Figure 10. Comparison of the self-similar shape of the splash (i.e., with all length scales rescaled 
by r 2//3 ) for four different values of the Weber number. Differences in the rescaled shapes are 
only observed for the lower values of the Weber number, see for example the shape for We = 68 
and We = 9 in figure [5] For high Weber numbers ( We > 100) the rescaled shape hardly changes 
as a function of We. 
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Figure 11. Length scales and velocities in the splash rescaled by r 2//3 and t -1//3 respectively. 
After proper rescaling, all the plotted values become independent of time, which indicates the 
self-similar behavior. The splash width wp defined in figure [5] position and vertical velocity at 
the point P, up, as a function of time. As expected from the analysis, the horizontal position 
x p of point P as well as the width w p of the splash at point P are proportional to r 2//3 . Due 
to the boundary condition at the disc, a constant B u is added to the vertical velocity u p and 
a term B z times r is added to the vertical position z p of point P (see text). All data are for a 
simulation at We — 274. 
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Figure 12. The positions x p and z p , the splash width w p , and the vertical velocity u p , for 
different values of the Weber number. Clearly, all the rescaled values are independent of the 
Weber number for We > 100, consistent with the fact that the self-similar shape of the splash 
becomes independent of We in this regime (see also figure [10)) . 



closer look at the breakup of the splash into a crown. With the crown breakup, we are 
referring to the ejection of droplets from the tip of the splash. 

In order to determine when the deceleration is large enough to provoke the growth of 
perturbations, we define a local Bond number at the tip of the splash using the radius of 
curvature and the acceleration felt by fluid particles in this region of the wave, namely, 
Atip — g (see equation (jl.ip ) , where 



-4 



tip 



d Utip 
dT 



(4.1) 
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Figure 13. Influence of the node density on the tip of the splash, for We = 274. (a) The 
radius of curvature Rc at the tip has a minimum value comparable to the minimum distance 
between the nodes. The radius of curvature increases as the cylindrical rim grows in size, until 
the rim pinches off, and the radius of curvature starts again at its minimum value. This process 
is sensitive to the node density, but has no influence on the solution away from the tip, and is 
not present when We is small enough, (b) The cylindrical rim at its maximum size before it 
pinches off, for a minimum node distance of 0.001. (c) The same simulation as (6), just after the 
rim has pinched off. The rim is removed from the simulation directly after the pinch-off because 
it does not represent the physical situation (see text). 



with —dUup/dT > 0. However, the determination of the local Bond number at the tip 
of the splash wave is not an easy task. Indeed, from the numerical simulations Bou p can 
be obtained very accurately, but only in an indirect manner for most Weber numbers. 
The reason for this is that the tip of the splash is unstable in the numerical simulations 
for We > 30, as is explained in more detail later in this section. Alternatively, the local 
Bond number at the tip can be determined directly from the experiments, but with a 
relatively large uncertainty. Both the experimental and numerical method of determining 
the local Bond number will be explained now. 

We obtain experimental values for Bo t i P by tracing the tip of the splash, within a 
time interval Ti with a duration of typically 3 ms, in high-speed movies taken at 5400 
frames per second. A second-order fit to the position data versus time gives us the mean 
deceleration of the tip A t i p . We determine the radius of the rim Rc by measuring it 
graphically at the beginning and at the end of Ti, giving us a minimum and maximum 
value for Rc within the time interval. The experimental values for Bou p are shown in 
figurc[l4]in black dots. The error bars are obtained by using the minimum and maximum 
of Rc, because the radius of curvature is the dominant source of uncertainty, due to the 
squared appearance of Rc in the definition of the Bond number (|1.1[) . 

In the boundary integral simulations, the radius of curvature and deceleration can 
be determined accurately at any moment of time. The tip of the splash, which is the 
region in which we are interested, is however only stable in a limited range of low Weber 
numbers. For higher Weber numbers ( We > 30) a neck is formed below the rim, leading 
to its pinch-off from the rest of the splash wave (see figure n"37 6-c)1fl. As soon as the 
pinch-off has occurred, a new neck forms which pinches off a bit later. This induces 
a series of pinch-offs in the numerics, which are unphysical for a number of reasons. 
The first is that, due to axisymmetry, the droplet is tubular, much unlike the droplets 

f This pinched-off donut-shaped volume of fluid is subsequently removed from the simulation 
to prevent it from crossing the free surface in some point which would cause the code to crash. 
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that are generated in experiments. The second is that the details of the pinch-off are 
strongly dependent on the node density that is used in the BI simulations, as is shown 
in figure 1131 and therefore needs to be qualified as a numerical artifact. The series of 
pinches influence both the length scale and the deceleration of the tip, such that it is 
not possible to determine the local Bond number at the tip of the splash for high Weber 
numbers in the simulations. Although the tip of the splash is unstable, we can show that 
the rest of the splash is not influenced by these numerical artifacts. For this reason, we 
determine the local Bond number at point P, Bop, using the width W p (see figure [5]) as 
the length scale and calculate the deceleration at t he same point, namely, 

Bo P =P(-^- 9 )wl (4.2) 



cr V dT 

In figure [T4l the different methods of measuring the local Bond number are put together. 
Using the lower reference point P on the splash shows that for We > 100, the local Bond 
number is proportional to We (blue squares, figure H4)). and we find that it is constant 
in time. Comparing Bop with the experimental values of Bou p however, shows that 
the local Bond number is overestimated if we use Bo p. But calculating the local Bond 
number in the numerical simulations at different positions on the splash wave results 
only in a vertical shift of the points in figure HU Indeed, figure shows that the value 
of the local Bond number evaluated at the point Q on the splash wave where the slope 
is -2 (thus, point Q is closer to the edge of the rim than point P), is proportional to Bop 
defined in equation (|4.2j) . This result is a consequence of the fact shown in figure QT] that 
the dependence of the lengths and velocities in the splash wave are not dependent on 
the Weber number for sufficiently large values of this control parameter. This evidence 
indicates that we can calculate the local Bond number at different positions by simply 
multiplying Bop by a constant. For Bou p , we determine this constant to be ~ 0.1 using 
the experimentally measured local Bond number at the tip of the splash. Also notice 
in figure [14] that the value for the constant ~ 0.1 is also consistent with the values of 
Boup numerically calculated for low values of the impact Weber number. The local Bond 
number at the tip that we have deduced from Bop is shown in figure [T4l with red circles. 
Clearly, the transition to breakup into a crown occurs when the local Bond number at 
the tip of the splash is of order unity. 

The proportionality of the local Bond number with the Weber number and its indepen- 
dence of time is not unexpected, because it results from the self-similarity solution that 
becomes independent of We for large Weber numbers, as shown in figure Q21 Using the 
self-similar scaling for distances, velocities and accelerations, we write the deceleration 
and typical length scale in dimensional form as: 

A UP = ^-C r-*/3 (4.3) 
tip, 

and 

Be = RdCrt 2 / 3 , (4.4) 

where C a and Cr are dimcnsionless constants, independent of time and Weber number 
for We > 1. Substituting (|4.3|) and (|4.4j) in (|1.1|) . and using the fact that A tlp > g for 
large values of We gives: 

Bo Up * C a C 2 B pV ° RD = C a C 2 p We , (4.5) 
a 

clearly showing that the local Bond number is independent of time and proportional to 
We. 
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Figure 14. The local Bond number measured in different ways as a function of We. The blue 
squares are the local Bond numbers determined in the simulations at point P (Bop), the black 
triangles correspond to the local Bond number at the point Q {Boq) of the inset and red circles 
are these same values of Bop multiplied by 0.1. This multiplication factor has been found by 
matching to the experimental values of the local Bond number at the tip (black dots). The 
green diamonds show the local Bond numbers measured at the tip of the splash in the BI 
simulations, which is only possible for We < 30 (see text). The shaded area indicates the range 
of experimental conditions in which we observe the crown breakup transition. Note that this 
transition occurs when the local Bond number is of order unity (dashed horizontal line). 



Using the proportionality Bo t i P oc We and the crown breakup condition that Boti P is 
of order unity around the transition, we can define a condition for the crown breakup 
transition based on We that does not involve the local Bond number. Such a condition is 
preferable, because whereas the local Bond number is difficult to determine, the Weber 
number is directly given by the experimental conditions. There indeed exists one critical 
Weber number We cr i t above which we always observe breakup of the splash into a crown, 
as can be seen in table Q] The table shows the values of We, Fr and Re at the crown 
breakup transition^]. Note that only the Weber number always has approximately the 
same value at the transition, showing that the Weber number is indeed the relevant 
parameter to indicate the transition to crown breakup and droplet ejection. 



| Clearly, We, FY and Re can be related after calculating the disc velocity Vd using the disc 
radius Rd and either of the three dimensionless numbers in table [1] 
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Rd (mm) We cr it Fr Re 

15 145 ±5 4.77 1.26 ■ 10 4 

20 135 ± 19 2.51 1.40 ■ 10 4 

25 134 ± 11 1.59 1.56 ■ 10 4 

30 142 ±19 1.18 1.77 • 10 4 

Table 1. Transition to crown breakup of the splash for different disc radii. The value of the 
critical Weber number does not differ appreciably with the size of the disc, where the Froude 
number and the Reynolds number at the transition show a clear dependence on Rd- 



5. Air cushion effect: force on the disc and drop size 

The previous theoretical derivation of the self-similar solution predicts decelerations 
within the splash wave proportional to oc i~ 4 / 3 , with t the dimcnsionless time after 
impact. The huge decelerations expected from the analysis for t <C 1 arc not observed 
in a real experiment since there are other physical effects, not taken into account in 
the previous theoretical description, that break the self similar structure of the flow 
field during a short time interval after impact. This is clearly seen in figure [3l where 
it is depicted that the splash wave initially propagates radially outwards at a velocity 
constant in time, and not at a diverging velocity proportional to oc t -1 / 3 , as is predicted 
by the self-similar solution. The origin of the cutoff time or, equivalently, the origin of the 
cutoff length below which the self similar solution is no longer valid, could in principle 
lie in any of the following three sources: (i) the fact that the disc edge is rounded, (ii) 
the effect of liquid viscosity, and (iii) the air layer between the impacting disk and the 
liquid. 

First we assess the fact that the disc edge is not perfectly sharp, but possesses a finite 
curvature. Since, in view of figure [5] the self-similar solution is only valid for t > to exp ~ 

2 /3 

0.03, we can estimate the cut-off length X Q = xqRd ~ RDt 0exp ~ 1.9 mm. This is 
comparable with the disc thickness and clearly much larger than the length associated 
with the finite radius of curvature at the disc edge which we estimate to be ps 50 /Ltm. 
To evaluate the effect of liquid viscosity we estimate the viscous cutoff time, tov, whose 
magnitude is determined by equating the boundary layer thickness, S, to the width of 
the splash region at the instant of time when the self similar solution begins to be valid, 
namely, 



hxRptov 2/3 .1/6 D -1/2 D _ 3 1n -l0 (r n \ 

° ~ \ 77 — ~ tin *o« ~^ Vu ~ & e o tov ~ Re D ~ 10 . (5.1) 

V P V D 

Clearly, this value is much smaller than to exp and viscosity can therefore not account for 
the observed cut-off time. 

The remaining physical effect is the finiteness of air de nsity, which is known t o be re- 



spons ible for the so called air cushion effect, analysed bv lHowison et al\ (|199lf k I Wilson 



( 19911 ) for the case of two dimensional geometries. Figure [T5k illustrates the disc ap- 
proaching the free interface forcing the air between disc and liquid surface to flow radi- 
ally outwards. Using continuity and neglecting compressibility effects, the radial velocity 
field of the gas, V g (R,T), can be expressed as a function of (dimensional) time T and 
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Figure 15. a) The downward motion of the disc forces the gas to flow radially outwards. The 
overpressure needed to accelerate the gas deforms the interface, b) The linearity of the Laplace 
equation permits to express the liquid potential at R = for a generic value of dH/dT as 
= 0, T) — —aRo dH/dT, with oRdVd the value of the potential at R = obtained from 
the solution of the Laplace equation subjected to the boundary conditions depicted in this figure. 
Here, d^/dn indicates the normal velocity at the interface. 



the vertical deformation H g (T) of the free surface below the disc as 
-irR 2 [Vd - ] + 2ttR [H g + V D (T s - T)} V g =0 



V„ = 



dT . 

R {V D -dH g /dT) 



(5.2) 



9 2 H g + V D (T S -T) 

where we assume that H g in the region below the disc does not depend on R. In 
T s is a constant which will be arbitrarily fixed to T s — Rd/Vd since its precise value 
possesses a negligible influence on H g (T). Using the expression for the flow field given 
in (|5.2p . the gas pressure dependence on R and T can be determined integrating the 
momentum equation 

"§ + ^T + < 5 ' 3 > 

which yields the following expression for P g (R,T): 

3 2 p2] f Vd- dHg/dT \ 2 Rj-R 2 1 d* H g 

^9 ^ s P 9 {n D n ){ Ha + y D{Ts _ T) J Pg 4 H g + Vd(T s —T) dT 2 

(5.4) 

The desired equation for H g (T) is determined from the Euler-Bcrnoulli equation par- 
ticularized at R — 0, namely, 



<9$ , ,1 ( dH^ 



P^ R =^ T ) + lp\^i) +P g (R = 0) = P a , (5.5) 

where P a indicates the atmospheric pressure. Indeed, assuming that the normal velocity in 
the liquid below the disc does not depend on R, the potential $ at R = can be expressed 
as a function of H g since, by virtue of the linearity of the Laplace equation, $(i? = 0) = 
—aRndHg/dT, with a = 2/n the value of the potential at R = corresponding to 
the solution of V 2 $ = subjected to the boundary conditions sketched in figure ll~5b 
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(see Iafrati k, Korobkir] ( 2011 )). Finally, the substitution of equation (|5.4j) into (|5.5 



and subsequent non-dimcnsionalization provides the following differential equation for 
h = Hg/Ri)-. 

lp g 1 \<Ph ifdhV 3p g f l-dh/dt \ 2 = 
4 p t s -t + h) dt 2 2\dt J 8 p \h + t s -t J ' { ' 

with t = TVd/Rd and similarly for t s . 

Let us point out that the above analysis rests on the assumption that viscous effects 
in the gas are unimportant, and this is indeed the case since 

Pg VrgHg _ PgVpR 
Pg 2p g 



Re g = Pg V ; 9 9 = ^lE^L „ O(10 3 ) , (5.7) 



where use of the continuity equation f|5.2[) has been made. The estimation in equation 
(|5.7[) . which suggests to neglect viscous effects as a first approximation, is valid while 
Vd — dHg/dT ~ 0(Vd), a condition which is fulfilled even when the disc position is 
below the unperturbed free interface, as it will be shown below. 

The integration of equation (|5.6|) . subjected to appropriate initial conditions, provides 
with the time evolution for h, dh/dt, the dimensionless gas pressure p(r = 0) = P g /(pVp) 
and the height of the gas layer hi = t s — t + h shown in figures fTBTITTl It is depicted in 
figure fTBTa) that the interface is only appreciably accelerated for t ~ 1 i.e., when the 
disc is sufficiently close to the undisturbed interface. Moreover, only for times such that 
t > tend with tend = 1.01, namely, when the disc position is only slightly below the 
level of the undisturbed interface, dh/dt ~ 1 and h\ ~ 0. Therefore, the liquid normal 
velocities below the disc are smaller than the disc impact velocity while t < t en d and 
equal to the disc velocity for t > tend- This fact implies that the speed at which the 
splash wave propagates must be, during a time interval t — t s < 0.01, smaller than the 
one predicted by the previous potential flow numerical simulations where the effect of air 
was neglected. This is just what figure |4] shows: the splash wave calculated numerically 
needs 0.6 ms less than the real wave to reach the wave radial position at T = 0.2ms 
(r = 0.1); from that spatial position onwards, the wave speed calculated numerically is 
identical to that measured experimentally. We can conclude that the air layer below the 
disc breaks the self similar structure of the splash wave during t en d — t s = to ga s — 0.01, 
a time interval very similar to the magnitude of the virtual origin in times, to exp = 0.03, 
deduced in figure [BJ 

To check the validity of the model presented above, we have carrie d out two-fluid 



poten ti al flow numerical s i mulat ions using the numerical code detailed in I Gordillo et al 



( 2007t ): iGekle fc Gordillo! ([201 lh . Both the density ratio p g /p and the impact Weber 



number have been varied in wide ranges of values. Figure [16] shows that the time evolution 
of the dimensionless height h is almost perfectly reproduced by our theory. The shape 
of the entrapped bubble calculated numerically for different values of the density ratio, 
depicted in figure [TTl shows that the interface accelerates beneath the disc uniformly 
downwards, except in a very localized region near the disc edge, where the interface 
accelerates upwards as a consequence of Bernoulli suction (or, equivalently, due to the 
growth of a Kelvin- Hclmholtz instability), causing the the entrapment of an air pocket. It 
is also shown in figure [T7] that the maximum height of the gas layer calculated numerically 
compares very favorably with the one predicted theoretically in spite of the fact that our 
model, which was built under the assumption that h docs not depend on r, is unable 
to capture the upwards acceleration of the free interface. Let us point out here that the 
air cushion effect described above shares many similarities with the bubble entrapment 
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Figure 16. (a) Using continuous lines, we plot the time evolutions of the gas layer depth h, the 
liquid velocity at r = 0, dh/dt, gas pressure at r = 0, p g , and the distance from the disc bottom 
to the liquid interface, hi = to — t + h for p g /p = 1.2 x 1CF 3 and a — 2/ty. The inset shows the 
time evolution of the different variables for instants of time right after the disc bottom is below 
the undisturbed free interface. The discontinuous line represents the function h(t) obtained 
from the two-fluid potential flow code, which is almost perfectly reproduced by the integration 
of equation ()5.6|) . (b) Time evolution of the functions representing the dimensionless force on the 
disc fo (dashed lines) and mo (continuous lines) defined, respectively, in equations (|5.9p - (l5.10p . 
The insets show that due to the effect of air the force is not infinite any more, but reaches a 
well-defined maximum slightly after the disc bottom reaches the level of the unperturbed free 
interface and that the maximum value of the dimensionless change in momentum, defined in 
(pTOj) . is m D (t = 1.01) ~ 0.315 ~ O(l). 
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Figure 17. Continuous lines indicate the dependence of both h and dh/dt on the density ratio 
p g f p at t = 0.99 predicted by the integration of equation (|5.6[) . Diamonds indicate the numerical 
values of the air layer thickness at r = when the first contact between the liquid and the disc 
bottom occurs, which occurs near the disc edge, as shown by the different insets. Numerical 
simulations reveal that, except close to the disc edge, the air layer depth is mostly uniform. 





processes taking p l ace ei t her during the impac t of drops against solid o r liqui d surfaces 
Thoroddsen et all (120031); iMandre et all (120091) ; iDuchemin fc Josserandl (|201l[) or in the 
collapse of bubbles lGo"rdillo fc Fontelosl (|2007t ); iGordillol (|2008l) . It should be mentioned 
here that, while the gas flow is well described by the lubrication approximation in the 
impact of drops, viscous effects are negligible in our case. 
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Figure 18. Comparison between the time evolution of the horizontal position xp of point P for 
We — 274 in (i) a single-fluid potential-flow numerical simulation in which the liquid normal 
velocities beneath the disc is provided by our theoretical prediction (function dh/dt in figure 
116a . blue circles) and (ii) the experiments (green diamonds). Note that the experimental data is 
the bare data of figure|6]3 without the time shift. The conclusion is that the simulations based on 
our theoretical analysis excellently reproduces the experiments, with no adjustable constants. 

The critical role played by the thin layer of air located beneath the disc on the splash 
wave speed can be demonstrated in an even more convincing way. To this end, we have 
carried out single-fluid potential-flow simulations in which the effect of air is introduced 
by replacing the constant impact velocity boundary condition by the function dh/dt 
represented graphically in figure [Tor. Figure (fP%|) shows that this new type of numerical 
simulations, which indirectly retain the air cushion effect through its influence on the 
normal velocity beneath the disc, perfectly reproduces the radial propagation of the 
splash wave and, thus, avoids the need of using an artificial time shift to match simulations 
with experiments. 

Once the critical role played by air in the impact process has been shown to be very 
well described by our theory, we determine the force on the disc, Fjj(T), by integration 
of the pressure field given by equation (|5.4[) . to give 



F D (T)=pTrR 2 D V*f D (t), 



(5.8) 



where the function 




h + t s - t dt 2 



2 d 2 h 



(5.9) 
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is represented graphically in figure [HJb) for the particular case of p g /p = 1.2 x 10~ 3 . 
Figure 116( b) shows the remarkable result that fo (t) reaches a maximum when the disc 
bottom is below the level of the unperturbed free surface but before the disc bottom 
touches the li quid. Let us point ou t here t hat t he force on an impacting disc was already 



calculated by llafrati fc Korobkinl (|2008l 120111 ). who neglected the effect of air in the 



analysis, a starting hypothesis leading to the force on the disc to tend to infinite during 
the initial instants of the impact process. The fact that our theory predicts a maximum 
force on the disc, indicates that the air cushion effect described above softens the impact 
process, thus avoiding the diverging force deduced from the analysis when air is removed. 
Our theory also allows to determine the change in the momentum of a disc falling freely 
and impacting against a liquid free surface, 

r T 

AMfl(T) = / F D (T')dT' = P TrR 3 D V D m D (t), (5.10) 
Jo 

where the function rriD(t) defined in equation (|5.10[) is also represented graphically in 
figure ITBT b). Since ra,o(t = 1.01) ~ 0.315, the change in velocity A Vn of a free falling 
disc of density p s and thickness H impacting normally onto water-air interface {p g /p = 
1.2 x 10~ 3 ) can be estimated as: 

Ps ir R% H Ps AV D ~ 0.315tt pV D R 3 D -> ~ 0.315-^^ . (5.11) 

V D Ps H 

The result expressed by equation (|5.1ip indicates that, since for most solids 1 < p s /p < 
20, the huge deceleration experienced by the disc could reduce the impacting velocity 
an amount ~ 0(Vd) before it touches the interface in the case of Rd/H > 10. Let us 
finally point out that the validity of our study is restricted to those cases in which the 
disc impact velocity is sufficiently low so as to neglect gas compressibility effects. 

5.1. Drop size 

It was demonstrated above that, in the range of times in which the self-similar solution 
is valid (i.e., when tQ exp « ( « 1), the local Bond number is constant in time. Indeed, 
during this time interval, distances and decelerations in the tip region are proportional 
to Rc/Rd = rup oc (t — toexp) 2 ^ 3 and au p oc (t — t,Q exp ) ~ 4 / 3 respectively and, thus, 

p Rd 

Boi oca i = — S. a tlp r 2 u cx We. (5-12) 

a 

Therefore, whenever Boi oca i > 1 or equivalently, if We > We cr u, it is expected that 
when the rim decelerates i.e., at any instant of time such that t > to exp , the rim develops 
a R ayleigh- Taylor instab i lity w ith a initial length scale d < r t i P given by the condition 



(see IVillermaux fc Bossal (1201 lh ) 



pV D RDdmaxJipd , _i/ 2 -1/2 Cr 1 o\ 

'— ~ 1 -> d ~ We I a ma ' H , (5.13) 

(7 

where it has been taken into account that the maximum deceleration at the tip of the 
splash wave is such that a ma x,tip 3> gR/V^. Equation (|5.13|) expresses that the char- 
acteristic length scale of the fingers, dRu, is set by the condition that their weigh per 
unit length, namely, pV^Ro atipd 2 , equals the surface tension force per unit length, a. 

The maximum deceleration scales as a ma x,tip x toexp ~ ^Ogas anc ^' therefore, the initial 
wavelength of maximum growth rate is given by 

d^tH^Wer 1 / 2 , (5.14) 
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Figure 19. a) Snapshot showing the crown splash and how the edge finding routine identifies 
the ejected drops. Drop diameter is experimentally determined from the analysis of this type of 
images, b) This figure compares the diameter of the ejected drops with the theoretical prediction 
given by equation (|5.16[) for two different values of b. 

where use of equation (|5.13|) has been made. These initial perturbations give rise to the 
formation of fingers of increasing width and length which break due to the action of 
surface tension forces in a characteristic capillary time given by 

tca P = V D T cap /R D = V D /R D ( P R 3 D d 3 /a) 1/2 = We^d 3 / 2 cx We- 1 / 4 t 0exp , (5.15) 

where use of equation (|5.14[) has been made. Consequently, the breakup time will be 
given by tbreak oc (to exp + t cap ) and the dimcnsionless size of the drop formed, which is 
of the order of the width of the rim (see figure the experimental evidence in figure [T9k ). 
is given by 

ddrop ~ w rim cx t 2 b H ak cx tlH p (1 + t cap /t 0exp ) 2/3 cx tUlj, (l +We _1/4 ) 1 , (5.16) 

with b an adjustable constant. Figure 119b shows a reasonable agreement between the 
experimental drop size measured from the analysis of the images of the type depicted 
in figure Q] and the drop diameter predicted by equation (|5.16|) . Let us point out that 
figure [T9b reveals that the drop size is hardly dependent of the impact Weber number. 
This fact indicates that the characteristic length scale of the ejected drops is fixed by the 
air bubble entrapped underneath the disc, being this the physical mechanism underlaying 
the experimentally observed cut-off time, to exp . 

6. Conclusions 

We have studied the formation of a splash wave and the transition to the ejection of 
droplets after the impact of a disc on a liquid by using boundary integral simulations, 
experiments and theoretical analysis. Only by combining these three methods, we have 
been able to analyze the full problem. Although all the different approaches used possess 
limitations either in accuracy, stability or mathematical formulation, each of them is 
able to reveal specific parts of the problem unaccessible by the other methods. Using the 
overlapping parts however, we accurately demonstrated the validity of each of the three 
different approaches used in this study. 

By approximating the time just after disc impact (t <C 1) as a 2-dimcnsional potential 
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flow problem and neglecting the effect of air, we have shown that there exist self-similar 
solutions of the second kind for any value of the Weber number. These self-similar solu- 
tions exist because the matching of the incrtial terms in the unsteady Bernoulli equation 
to the far-field boundary condition for the potential gives the same scaling powers as the 
matching with the surface tension term. When the Weber number is increased to large 
values (> 100), the shape of the splash becomes independent of the Weber number. Both 
predictions are confirmed by boundary integral simulations, by rescaling the calculated 
shapes for a wide range of Weber numbers. We found the correct scaling for both the 
lengths (~ i 2 / 3 ), and velocities (~ i -1 / 3 ) at different positions in the splash. In the ex- 
periments, the same scaling is found after introducing a small time-shift to exp which we 
demonstrated is associated to the finiteness of gas density. Indeed, in spite of its large 
density contrast with water or the solid, the very thin air cushion entrapped between the 
impacting solid and the liquid reveals that gas critically affect the velocity of propagation 
of the splash wave as well as the force on the disc, Fn- Th e presence of air avoids th e 
singularity in Fd predicted when air effects are neglected (see llafrati fc Korobkinl (|201lh ) 
and our theory predicts that the maximum of Fn is reached before the disc touches the 
liquid. 

We have shown that the transition to droplet ejection (crown breakup transition) is 
caused by a Rayleigh- Taylor instability. We experimentally determined that the local 
Bond number based on the tip deceleration and the tip width, is of order unity at the 
splash transition. Theoretical analysis predict, and numerical simulations confirm that, 
Boup depends linearly on the Weber number for We > 100. From this linear dependence 
between Bou p and We, we have concluded that the splash transition can be identified by 
a critical Weber number, which we indeed find in the experiments. We have also found 
that the air cushion effect also sets the maximum deceleration experienced by the rim 
edge and have used this information to scale the sizes of the drops ejected for those exper- 
imental conditions in which We > We cr it — 140 and the impact Weber number is below 
the threshold for which surface sealing occurs. Finally, we have found the interesting re- 
sult that the thickness of the air layer entrapped beneath the disc sets the characteristic 
length scale of the drops ejected when the impact Weber number exceeds the critical one. 
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